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ABSTRACT 


A  class  of  numerical  dissipation  models  for  central-difference  schemes  constructed  with 
second-  and  fourth- difference  terms  is  considered.  The  notion  of  matrix  dissipation  asso¬ 
ciated  with  upwind  schemes  is  used  to  establish  improved  shock  capturing  capability  for 
these  models.  In  addition,  conditions  are  given  that  guarantee  that  such  dissipation  models 
produce  a  TVD  scheme.  Appropriate  switches  for  this  type  of  model  to  ensure  satisfaction 
of  the  TVD  property  are  presented.  Significant  improvements  in  the  accuracy  of  a  central- 
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I.  Introduction 


Central  difference  type  schemes  are  currently  being  applied  on  a  regular  basis  in 
the  solution  of  the  Euler  and  Navier-Stokes  equations.  A  numerical  dissipation  model  is 
included  in  these  schemes,  and  it  plays  a  crucial  role  in  the  determination  of  their  success. 
The  form  of  the  dissipation  model  is  quite  often  a  blending  of  second-difference  and  fourth- 
difference  dissipation  terms.  The  second-difference  terms  are  used  to  prevent  oscillations  at 
shock  waves,  while  the  fourth-difference  terms  are  important  for  stability  and  convergence 
to  a  steady  state.  There  is  a  constant  to  be  specified  for  each  contribution.  However, 
using  the  model  in  conjunction  with  appropriate  numerical  procedures,  these  constants 
can  usually  be  selected  and  maintained  for  a  fairly  wide  class  of  fluid  dynamics  problems. 
This  dissipation  model  allows  shock  waves  to  be  captured  with  smearing  over  three  to  four 
mesh  cells. 

Even  though  these  central-difference  schemes  have  proven  to  be  reasonably  effective 
in  many  cases,  there  are  still  strong  motivations  for  reducing  the  numerical  dissipation 
being  produced.  For  example,  by  appropriate  reduction  of  the  artificial  dissipation,  shock 
wave  representation  and  boundary-layer  definition  (especially  the  wall  shear  stresses)  can 
be  improved  on  coarse  meshes.  Such  improvements  in  accuracy  are  especially  beneficial 
for  complex  three-dimensional  flows,  which  can  demand  extensive  computational  effort. 
In  addition,  better  estimates  of  the  limit  of  infinitely  fine  mesh  values  of  aerodynamic 
coefficients  for  flows  with  shocks  can  be  obtained.  Also,  the  standard  model  has  difficulties 
in  hypersonic  flow.  Finally,  for  some  problems,  the  influence  of  numerical  dissipation  needs 
to  be  severely  limited  in  certain  smooth  regions  of  a  flow  field  (i.e.,  near  the  trailing  edge 
of  an  airfoil),  while  still  maintaining  stability  near  discontinuities.  This  difficulty  cannot 
generally  be  resolved  by  simply  reducing  the  global  constants  in  the  dissipation  model. 

One  can  appeal  to  ideas  from  upwind  schemes  to  improve  the  dissipation  model, 
especially  in  the  vicinity  of  shock  waves.  Upwind  algorithms  utilize  concepts  from  charac¬ 
teristic  theory  in  order  to  determine  the  direction  of  spatial  differencing.  They  have  been 
extended  to  systems  of  conservation  laws  using  such  approaches  as  the  flux-vector  splitting 
of  van  Leer  [1]  and  the  approximate  Riemann  solver  of  Roe  [2].  A  fundamental  feature 
of  these  schemes  is  the  use  of  a  matrix  evaluation  of  the  dissipation  either  in  the  implied 
or  direct  sense.  In  so  doing,  the  dissipative  terms  of  each  discrete  equation  are  scaled  by 
the  appropriate  eigenvalues  of  the  flux  Jacobian  matrices  of  the  Euler  equations,  rather 
than  the  same  eigenvalue  as  in  the  dissipation  model  employed  with  central-difference 
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schemes.  Also,  upwind  schemes  can  be  designed  to  have  the  total  variation  diminishing 
(TVD)  property,  which  prevents  the  occurrence  of  spurious  oscillations.  The  disadvantage 
of  these  schemes  is  that,  in  general,  they  increase  the  operational  count  for  processing 
mesh  points  by  about  a  factor  of  two  over  that  required  by  central-difference  schemes. 
One  would  certainly  like  to  more  closely  imitate  the  highly  desirable  behavior  of  the  up¬ 
wind  algorithms  near  flow  discontinuities,  and  at  the  same  time,  retain  the  more  efficient 
central-difference  scheme  over  significant  portions  of  a  flow  field.  In  addition,  one  would 
like  to  have  the  high  degree  of  numerical  efficiency  that  has  been  achieved  by  combining 
a  central-difference  scheme  with  a  Runge-Kutta  time  marching  algorithm,  which  includes 
residual  smoothing  and  multigrid  acceleration  techniques. 

The  primary  purpose  of  this  paper  is  to  construct  a  numerical  dissipation  model 
for  a  central-difference  scheme  that  has  both  the  properties  of  matrix  dissipation  and  of 
TVD.  As  a  starting  point,  we  consider  the  elements  of  a  widely  used  scalar  dissipation 
model.  Modifications  of  this  model  that  facilitate  accurate  viscous  flow  computations 
are  also  examined.  In  the  next  section  of  the  paper,  the  intimate  connection  between 
the  formulation  for  an  upwind  scheme  and  a  centered-difference  scheme  is  presented,  so 
as  to  establish  a  foundation  for  a  matrix  dissipation  model.  Then  a  theorem  is  proved 
that  provides  a  simple  sufficient  condition  to  determine  when  a  central-difference  scheme 
with  dissipation  terms  comprised  of  second  and  fourth  differences  is  TVD.  In  the  following 
section,  appropriate  flux  limiter  functions  consistent  with  the  central-difference  dissipation 
model  are  discussed.  A  multistage  time-stepping  scheme  used  in  applications  is  next  briefly 
described.  Finally,  numerical  results  are  shown  to  demonstrate  the  benefits  of  using  the 
matrix  dissipation  model.  Both  inviscid  and  viscous  transonic  airfoil  flows  are  computed. 

II.  Scalar  Dissipation  Model 

The  basic  elements  of  the  scalar  dissipation  model  considered  in  this  paper  were  first 
introduced  by  Jameson,  Schmidt,  and  Turkel  [3]  in  conjunction  with  Runge-Kutta  explicit 
schemes.  This  dissipation  model  has  been  used  by  many  investigators  [4-6]  to  numerically 
solve  the  Euler  equations  for  a  wide  range  of  fluid  dynamic  applications.  The  same  type  of 
dissipation  model  has  been  applied  to  alternating  direction  implicit  (ADI)  schemes  [7]  and 
LU  factored  implicit  schemes  [8].  Several  modifications  of  the  model  have  been  investigated 
in  [9]  and  [10]  in  order  to  improve  it  and  make  it  suitable  for  obtaining  accurate  and  efficient 
solutions  of  the  Navier-Stokes  equations.  In  this  section,  the  basic  model  and  important 
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modifications  are  briefly  reviewed. 

Consider  the  Euler  equations  in  the  form 

W*  +  fx  +  ffy  =  0,  (2-1) 

where  W  is  the  four-component  vector  of  conserved  variables,  and  /,  g  are  the  flux  vec¬ 
tors.  The  independent  variables  are  time  t  and  Cartesian  coordinates  (x,y).  If  (2.1)  is 
transformed  to  arbitrary  curvilinear  coordinates  £  =  £(x,y)  and  t)  =  y(i,y),  then  we 
obtain 

{J~1W)t  +  Fi  +  Gt,=  0,  (2.2) 

where  J-1  is  the  inverse  transformation  Jacobian,  and 

F=fyr,~gx  r>,  G  =  gxi-fyi. 

In  a  cell-centered,  finite-volume  method,  (2.2)  is  integrated  over  an  elemental  volume  in 
the  discretized  computational  domain,  and  J~l  is  identified  a .s  the  volume  of  the  cell. 
Equation  (2.2)  can  also  be  written  as 

J~lWt  +  AWi  +  BWr,  =  0, 


where  A  and  B  are  the  flux  Jacobian  matrices  defined  by  A  =  dF/dW  and  B  =  dG jdW . 

To  advance  the  scheme  in  time  we  use  a  multistage  scheme.  A  typical  step  of  a 
Runge-Kutta  approximation  to  (2.2)  is 

WW  =  w(o)  -ak£L  +  DrjG^k~^  -  ,  (2.3) 

where  and  are  spatial  differencing  operators,  and  AD  represents  '  ae  artificial 
dissipation  terms.  The  dissipation  terms  are  a  blending  of  second  and  fourth  differences. 
That  is, 


where 


r 

l 

(2.4) 

/.  . 

(2.5) 

(2.6) 
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and  A^,  are  the  standard  forward  and  backward  difference  operators  respectively  as¬ 
sociated  with  the  £  direction.  The  variable  scaling  factor  A  was  originally  chosen  as 

A<+*,j  =  2  +  Mt+i.y  +  (A^«,y  +  (A*f)*+i,y_  .  (2-7) 

where  A^  and  Xn  are  proportional  to  the  largest  eigenvalues  of  the  matrices  A  and  B.  The 
scaled  spectral  radii  A^  and  An  are  given  by 

A*  =  WVr,  -  vxn\  A  cyjx*  +  y2, 

A„  =  |w*€  -  uyt  \  +  C\Jx^  +  yj, 

where  u  and  v  are  the  Cartesian  velocity  components,  and  c  is  the  speed  of  sound.  The 
coefficients  and  e^4)  are  adapted  to  the  flow  and  are  defined  as  follows: 

et(+±  -•  =  kW  ma Vij,  i/i+ij,»i+2,j),  (2-8) 

*  '  2 


Ui,j  = 


_  P»+i,y  2 pij  +  pi~i,j 
Pt+X,y  +  2p»',y  A  Pi-l, y 


E%,>  =  max  [°'  ('t<<)  -  £.<+)i.;)]  ’ 


(2.10) 


where  p  is  the  pressure,  and  the  quantities  K.W  and  /c^4^  are  constants  to  be  specified.  The 
operators  in  (2.4)  for  the  r?  direction  are  defined  in  a  similar  manner. 

The  second-difference  dissipation  term  is  nonlinear.  Its  purpose  is  to  introduce  an 
entropy-like  condition  and  to  suppress  oscillations  in  the  neighborhood  of  shocks.  This 
term  is  small  in  the  smooth  portion  of  the  flow  field.  The  fourth-difference  dissipation  term 
is  basically  linear  and  is  included  to  damp  high-frequency  modes  and  allow  the  scheme  to 
approach  a  steady  state.  Only  this  term  affects  the  linear  stability  of  the  scheme.  Near 
shocks  it  is  reduced  to  zero. 

The  isotropic  scaling  factor  of  the  original  dissipation  model  as  given  in  (2.7)  is  gen¬ 
erally  satisfactory  for  inviscid  flow  problems  when  typical  inviscid  flow  meshes  (i.e.,  cell 
aspect  ratio  0(1))  are  used.  The  factor  can  produce  too  much  numercial  dissipation  in  the 
cases  of  meshes  with  high  aspect  ratio  cells.  This  is  also  an  important  consideration  for 
high  Reynolds  number  viscous  flows,  where  a  mesh  providing  appropriate  spatial  resolution 
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can  have  cell  aspect  ratios  O(103).  In  [9]  and  [10]  this  difficulty  is  remedied  by  replacing 
the  factor  of  (2.7)  with  the  anisotropic  one 

=  2  +  Wi+u. 

where 

(*c),\y  =  <Mr)  (A«)t,;’ 

<t>i,j{r)  =  1  +  ri,y  0<f  <1, 

and  r  =  A„/A^.  In  the  normal  direction,  one  defines 

(A’?)t  ,j  =  ( A»?  )  t.J  * 

Alternatives  to  the  switching  function  presented  in  (2.9)  have  been  investigated.  Cau¬ 
tion  must  be  exercised  in  the  selection  of  a  switching  variable.  If  a  quantity  with  the  same 
functional  dependence  as  entropy  (i.e. ,  p/p”1)  is  used,  sharper  shocks  can  be  captured  in 
viscous  transonic  flows.  However,  such  a  choice  can  result  in  a  loss  of  accuracy  for  the 
surface  shear  stress,  due  to  the  significant  variation  in  the  entropy-type  variable  across  the 
boundary  layer.  This  difficulty  can  be  removed  by  simply  multiplying  the  scaling  factor 
by  a  function  of  the  local  Mach  number  of  the  flow.  An  acceptable  modifying  function  has 
proven  to  be  (Ml/Moo)^,  for  some  (3  >  1,  where  Ml  is  the  local  Mach  number,  and  Mref 
is  a  reference  Mach  number  (i.e.,  free-stream  value  for  external  flows).  It  is  important 
to  note  that  the  entropy-type  function  is  generally  not  satisfactory  for  inviscid  flows.  In 
addition,  one  can  consider  the  sum  of  two  switches,  one  depending  on  pressure  and  the 
other  on  temperature  so  that  all  thermodynamic  changes  are  taken  into  account.  When 
introducing  the  matrix-valued  dissipation  it  will  be  possible  to  use  separate  switches  for 
different  characteristic  variables. 

The  treatment  of  the  artificial  dissipation  must  be  modified  at  the  boundaries  of  the 
physical  domain.  In  the  case  of  the  fourth-difference  dissipation  the  standard  five  point 
difference  stencil  must  be  replaced  at  the  first  two  interior  mesh  cells.  This  means  that 
one-sided  or  one-sided  biased  stencils  are  used  at  these  cells.  The  dissipative  character 
of  the  artificial  terms  is  important  because  it  influences  both  stability  and  accuracy.  For 
example,  if  the  dissipation  is  too  large  at  a  solid  boundary,  an  artificial  boundary  layer  is 
created  in  an  inviscid  flow,  and  the  effective  Reynolds  number  for  a  viscous  flow  is  altered. 
To  improve  accuracy  at  the  wall  boundaries  of  viscous  flows,  where  gradients  are  steep 
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due  to  physical  boundary  layers,  the  usual  fourth-difference  stencils  are  changed  in  this 
dissipation  model. 

Let  the  total  dissipation  for  a  mesh  cell,  in  the  direction  represented  by  the  index  j, 
be  denoted  by  dj.  For  simplicity  assume  that  Ae^4)  =  1.  Then, 

dj  —  dfj+^  ~  dfj 

where  the  dissipative  flux 

df,+i  =  (AHOi+j  -  2(A W)j+i  +  (A W)^ 

and  thus 

dj  =  (AlV);  +  }  -  3(A»')iti+3(AH');4  -  (AlV)^j  (2.11) 

with  the  index  i  for  AW  suppressed  for  convenience.  Consider  the  first  two  interior  cells 
adjacent  to  a  solid  boundary,  as  depicted  in  Figure  1.  If 


(AtV)i  =  (AW),  =(AW)k 


(2.12) 


then  (2.11)  gives 


d2  =  W4  -  2 W3  +  W2, 


(2.13) 


d,  =  W5  -  4Wa  +  5 W3  -  2W2. 


(2.14) 


These  boundary  stencils  are  fairly  standard  ones,  and  they  result  in  a  nonpositive  definite 
dissipation  matrix  for  the  system  of  difference  equations  [7].  An  alternative  form,  which  has 
reduced  the  sensitivity  to  solid  surface  normal  mesh  spacing  for  turbulent  flow  calculations 
without  compromising  stability  or  convergence,  is  given  by 


and 


(Aiy) i  =  2(A1F)|  -  (A1F)| 

(2.15) 

d2  =  W4  -  3WS  +  3VF2  -  Wu 

(2.16) 

d3  =  W5-  4W4  +  6W3  -  4W2  +  Wx. 

(2.17) 

This  boundary  condition  is  advantageous  if  the  mesh  is  fine  enough  to  adequately  represent 
the  laminar  sublayer  region  of  the  boundary  layer  (i.e.,  at  least  two  mesh  points  are  inside 
the  sublayer).  For  coarse  meshes  this  treatment  can  be  less  accurate  than  the  zeroth-order 
extrapolation  of  (2.12). 
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III.  The  Upwind  Connection 

Upwind  schemes  for  solving  hyperbolic  systems  of  conservation  laws  (i.e.,  Euler  equa¬ 
tions  of  gas  dynamics)  rely  upon  characteristic  theory  to  determine  the  direction  of  prop¬ 
agation  of  information,  and  thus,  the  direction  required  for  one-sided  differencing  approxi¬ 
mations  of  the  spatial  derivatives.  With  such  schemes  shock  waves  can  be  captured  without 
oscillations.  Thus,  a  successful  artificial  dissipation  model  for  a  central-difference  scheme 
should  imitate  an  upwind  scheme  in  the  neighborhood  of  shocks.  We  now  review  the 
connection  between  these  two  types  of  schemes. 

Consider  the  one-dimensional  scalar  wave  equation 


ut  +  aux  —  0 


un+1  =  uj  —  a- 
1  1  Ax 


(3.1) 


with  a  constant.  The  first-order  upwind  scheme  can  be  written  as 

|  Uj+ 1  -  Uj,  a  <  0 

Uj  —  Uy_  1  ,  a  >  0, 

where  all  discrete  quantities  are  evaluated  at  time  level  nAt  unless  otherwise  denoted.  The 
scheme  of  (3.1)  can  be  rewritten  as 


u 


n+ 1 


At 


At 


uj  ~  -  u/-i)  +  MiTxrK'+i  ~  2uj  + 


(3.2) 


3  J  2AxK  J  1  '2Ax 

which  now  ontains  a  central-difference  term  and  a  second-difference  dissipation  term.  Now 


consider  the  system 


ut  +  Aux  —  0, 


(3.3) 


where  u  is  a  N-component  vector.  The  system  case  can  be  converted  to  a  scalar  one  by 
diagonalizing  the  N  x  N  matrix  A  with  a  similarity  transformation 

A  =  AT, 

where  the  columns  of  T  are  the  right  eigenvectors  of  A.  After  diagonalizing  (3.3)  and 
applying  the  scheme  of  (3.2),  the  first-order  upwind  scheme  is  given  by 


tt"+1  =  Uj  ~  ^9AT^Uj  +  1  ~  U3~l)  +  ~  2uJ  +  UJ- l)’ 


At 


At 


2Ax 


2Ax 


(3.4) 


where 


A\  =  T\A\T~\  A  =  Diag  [|Ai|---|A 


N  I 
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The  generalization  to  a  system  of  conservation  laws  is  straightforward;  namely, 

ut  +  fx  =  0 

with  /  being  a  N-component  flux  vector,  and 

u]+l  0  +  ^  [|>4i+*|(ttJ+1  -uj)~  |^/-*|  (uj  ~  >  (3-5) 

where  the  Jacobian  matrix  A  =  df /du,  and  |A|  is  defined  as  for  (3.4).  The  matrix 
A,-:  i  can  be  computed  as  an  arithmetic  average  or  a  Roe  average.  For  transonic  steady 
flows  the  differences  are  negligible;  therefore,  we  use  the  simpler  arithmetic  average.  For 
hypersonic  flows  Yee  [11]  found  that  the  Roe  average  yields  better  results.  For  time- 
dependent  problems  the  Roe  average  also  seems  to  give  slightly  better  results. 

IV.  Matrix  Dissipation  Model 

We  now  extend  the  scheme  given  in  (3.5)  to  the  two-dimensional  equations  of  fluid 
dynamics.  In  particular,  consider  the  transformed  Euler  equations  of  (2.2)  with  the  Runge- 
Kutta  scheme  of  (2.3).  The  necessary  modification  to  the  contributions  for  the  £  direction 
of  the  artificial  dissipation  term  defined  by  (2.4)  is  to  substitute  \A\  for  the  eigenvalue 
scaling  factor,  A,  in  (2.5)  and  (2.6).  For  the  r?  direction,  £  and  |A|  are  replaced  by  77  and 
|f?|,  respectively.  We  next  define  explicitly  the  form  for  the  matrix  |A|.  Let 

A  =  Diag  [Aj  A2  A3  A3] 
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0  0  0  O' 

—ayq  a2  aya  2  0 

—a2q  H2(ly  a2  0  1 

—q2  q<n  qa2  0. 

—q  a  1  a2  O' 

— uq  ua  1  ua2  0 

—vq  va  1  va2  0  ’ 

—Hq  Hay  Ha2  0. 

0  0  0  0 
ay(j)  — CLyU  —  ayV  (ly 
a2<f>  —a2u  —a2v  a2 
q<p  —qu  —qv  q 

H  is  the  total  enthalpy,  and  4>  =  ( u 2  +  v2)/ 2.  Because  of  the  special  form  of  \A\  (i.e.,  each 
row  of  Ej  is  a  scalar  times  the  first  row)  for  any  A1(  A2,  and  A3,  an  arbitrary  vector  x  can 
be  multiplied  by  |A|  very  quickly.  That  is,  we  calculate  AJ+  1  ( u]  +  y  —  Uj )  directly  (see 

[12])  rather  than  calculate  AJ  +  i  and  multiply  a  matrix  times  a  vector.  The  matrix  \B\ 
is  computed  in  the  same  way  as  |A|  by  simply  replacing  £  with  r). 

In  practice  one  cannot  choose  A!,A2,A3  as  given  above.  Near  stagnation  points  A3 
approaches  zero  while  near  sonic  lines  Aj  or  A2  approach  zero.  A  zero  artificial  viscosity 
would  create  numerical  difficulties.  Hence,  we  limit  these  values  as 

|Ai|  =  max(|A1|,Vnp(A)),  p{A)  =  |  q\  +  c\Ja\  +  a\, 

|A2j  =  max(|A2|,  Vnp(A)),  |A3|  =  max(|A3|,  Vtp{A)), 

where  the  linear  eigenvalue  A3  can  be  limited  differently  than  the  nonlinear  eigenvalues. 
The  parameters  Vn  and  Vi  have  been  determined  numerically.  Various  values  have  been 
evaluated  by  comparing  their  corresponding  computed  solutions  on  the  basis  of  the  follow¬ 
ing:  1)  Sharpness  of  shock  waves  captured  (without  producing  oscillations),  2)  Convergence 
rate  of  numerical  scheme.  A  good  choice  for  Vn  and  Vt  is  between  0.2  and  0.3. 

We  have  thus  far  replaced  A,+  i. ;  in  (2.5)  and  (2.6)  by  a  matrix  while  leaving  the 
limiters  and  as  scalars.  One  can  also  introduce  and  into  the  diagonal 
matrix  A.  This  allows  different  limiters  to  be  chosen  for  different  characteristic  variables. 
For  example,  the  limiter  may  be  based  on  pressure  for  the  nonlinear  waves.  However,  the 
pressure  is  smooth  through  a  contact  discontinuity.  Hence,  a  switch  based  on  temperature 
may  be  more  appropriate  for  the  linear  wave.  One  could  also  use  different  mesh  scalings, 
4>{r),  for  the  linear  and  nonlinear  waves.  Also,  a  smoother  cutoff  [13]  can  be  introduced. 
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V.  The  CVD  Property 


Consider  the  one-dimensior>al  scalar  conservation  law 


^[“(*.01  +  J^[/(u(M))]  =  °> 


—  oo  <  x  <  oo,  t  >  0. 


Le„  u(t)  =  be  the  approximate  solution  of  (5.1)  and  consider  the  semidiscrete 

equation 


itvA *> +  ~  />-1  *  2L  h+t  Av^i  ~  Q>-iAvi-i. 

"  aT  n+^3,,;+i  ~ 


Auy+.  =  (Aw)y+.  =  vy+1(<)  -  vy(0; 


A3  is  a  third-difference  operatoi  defined  as 


A  Wy+i.  =  (A  v)y+i  =vj+2{t)  -3vj+1{t)  +3 Vj(t)  -wy_!(0. 

The  terms  on  the  right  hand  si^v.  of  (5.2)  represent  second-  and  fourm-difference  numerical 
dissipation  terms,  with  a  constant.  Define 


S;  +  i  =  sgn 


where  sgn  represents  the  signum  function.  Then,  a  modification  of  the  theorem  of  Tadmor 
[14]  is  given  by  the  following. 

Theorem:  The  semidiscrete  scheme  of  (5.2)  is  TVD  if 


A/y+1 

2~  sj-\  +si+|)j  Qj+k  -  5i+Jr(5j'+|  _  “J-  0'^v.  7 


Rj+  *  =  0  when  Sj+ 1  -  2sy+  *  +  «/--  1  /  0- 
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Proof:  Shifting  the  indices  by  one  in  (5.2)  and  subtracting  (5.2)  from  the  resulting  equation, 
we  obtain 


d  A  a.  1 

—  Av.  ,  i  + 


dt 3  +  *  2  Ax 

1 


A/j  + 1  AA'~? 


2Ax 


Ax 


Q,+ 3  Auy+3  -  2g;>,  A«y+*  +  Q;  - ■  . 

fii+|A\+|-2i2y+i  A3u.+  .  +^_.A3ur_._  - 


(5.3) 


Multiply  (5.3)  by  sj+lr  and  sum  over  all  j.  Note  that  sj+^  -  ±1,  so  s2j+^  1,  and 

s3  +  k  AvJ+k 


Av;  +  ±- 


We  then  obtain 


A 

dt 


X/  lAt;j  +  yl  ”  2AzXsf+J  p-J  sj+h  AvJ  + 1 


2  'Avi+k 


+ 


2^T,sJ+k  (SJ>§  -  2*J  +  *  +  SJ-0  +  * 

j 

Ov+I  “25>+l-  +  5>-2-)  ^+^3w>+y* 


Au/+i 


(5.4) 


We  stress  that  the  last  term  in  (5.4)  will  not  help  for  TVD.  Its  purpose  is  to  eliminate  high 
frequencies  and  accelerate  convergence  to  a  steady  state.  Hence,  we  want  this  contribution 
to  be  zero.  This  can  be  accomplished  if  we  demand  either 


sj+§  ~2sJ  +  ±  +5 


0  2 


or 


*y+*=0’ 


(5.5) 


(5.6) 


i.e.,  condition  (b).  We  are  left  with 


3;  (TV)  =-5sE|-W*<+ 1 

3 


AA 


-H 


Av3+k 


-  sj+^s3+i  ~  2s/+  k  + 


+  5 


Av,- 


i  +  T 


(5.7) 


where  TV  denotes  the  total  variation  as  given  by 


Thus,  a  sufficient  but  not  necessary  condition  that  the  total  variation  not  increase  is  that 
the  term  of  (5.7)  in  brackets  must  be  positive.  This  means  that 


{3]  +  §■  2sj  +  k  +  s:-k)  Qj+k-  si+h  +  !  si~k)  At/  \ 


j+i 


(5.8) 


Since  s2  ,  =  1,  (5.8)  is  equivalent  to  condition  (a)  of  the  theorem.  Defining 

3  '•  2 


X]  =  1  - 


3- 


1  Si4- 

2f  •?  + 


(5.8)  can  be  written  in  the  form  given  by  Tadmor;  namely, 

A/+i. 

[X]  +  Xj+i)Q]+^  >  (xj  -  xy+i)  . 

Several  remarks  concerning  this  theorem  are  in  order. 

Remarks 

(i)  When  no  extrema  are  present  locally  (sy_i.  =  Sy+ 1  =  s;+.|),  both  conditions  of 
the  theorem  are  trivially  satisfied  for  all  QJ  + 1  and  Rj+i-  For  such  regions  we  want 
Qj  +  l  =  O(Ax)  for  second-order  accuracy  and  R}-+  ±  to  be  chosen  so  that  high  fre¬ 
quencies  are  damped. 

(ii)  If  Sj_  l.  =  Sj+  3  =  —  Sj+  l.  (a  local  oscillation  at  Xj+  i.),  we  require 


Qj+j  —  ^3  +  k  — 

(iii)  If  5y_i  =  Sy+i  =  —  -sJ+3  (a  local  extremum  at  xJ  + 1), 


* 


■A/; 


+  : 


5-  -  Av 


R;+  i.  -  0. 


(iv)  If  Sj_  i  =  —  s;-+  l.  =  —Sj+  3  (a  local  extremum  at  z;), 


Q;+i>  Rj+i  =  o. 

7  2  Av,-,  i 

■>  T  2 


It  follows  from  these  inequalities  that  <?J+i  can  be  negative.  As  far  as  total  variation  is 
concerned,  central  differences  are  not  nondissipative.  That  is,  they  can  either  increase  or 
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decrease  the  total  variation.  In  cases  where  central  differences  decrease  the  total  varia¬ 
tion,  Q ■ ,  i  can  be  negative.  For  systems  of  equations,  especially  in  multidirnensions,  this 
behavior  can  sometimes  lead  to  difficulties.  Hence,  we  demand  the  stronger  condition  that 


Qj+k  > 


Afj+k 


For  systems  this  condition  is  replaced  by 


> 


5 


where  A  =  df  /du ,  and  the  average  at  j  +  i  is  constructed  by  the  technique  of  Roe  [2]. 

VI.  Flux  Limiters 


In  this  section  switching  functions  are  introduced  that  force  the  scheme  to  automat¬ 
ically  satisfy  the  inequalities  presented  in  Section  V.  These  switches  are  required  to  be 
smooth  so  that  limit  cycles  are  not  experienced  when  marching  in  time  to  obtain  a  steady- 
state  solution.  The  discrete  switching  functions  are  defined  as 


Qj  +  l  = 


Ri+k  =  r  ' 


with  0  <  ip  <  1, 


(6.1) 


where  xp  ~  1  near  extrema,  so  that  the  inequalities  are  satisfied,  and  xp  —  0(Ax)  in  smooth 
regions  of  the  flow  field.  Conversely,  T  =  0  at  extrema,  and  T  =  1  in  the  smooth  regions. 
The  functions  t/>  and  T  of  (6.1)  can  be  defined  in  terms  of  a  limiter  function.  Let 

r  =  v3  Zlizl  - 


Vj+1  -  vj 


The  van  Leer  flux  limiter  is  given  by 


ife.  ^>° 


0,  r  <  0. 

From  Sweby  [14]  it  is  straightforward  to  see  that  for  any  flux  limiter 


V’y(r)  =  1  -  <P}(r). 
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Since  we  want  xpj(r)  to  be  positive,  this  relation  is  redefined  as 


0j(r)  =  |1  -  <pj{r)\. 


Then,  for  the  van  Leer  limiter, 


so  0  <  ip j  <  1  and  define 


1  -  r 

TTm 


^  r>0 


1,  r  <  0, 


0J  +  *.  =  max(0y,  ipj+\)- 


We  now  show  that  the  inequality  (5.8)  is  satisfied  with  this  xp.  By  the  first  remark  of 
Section  V,  the  inequality  is  satisfied  when  s-  i  =  s,  ,  t  =  s,-,  3,  and  so  in  this  case  we 

*  J  2  J  i  2  7  '  2 

only  need  Q  =  0(A:r).  Since  r  >  0  in  this  case, 

«'>  =  TT7- 

In  addition,  for  smooth  regions  of  the  flow  field, 

r  =  1  +  O(Ax), 


and  thus 


fy+±  =  O(Ai),  Qj+  .  =  O(Az). 


Next,  consider  the  case  s,_i  -  — s,  ,  j_.  This  implies  r,  <  0,  and  so,  xpj  =  1.  Moreover, 

J  2  •*  '  2  J  J 


S.L 
J  2 


sj+k  =>ip,  =  1  =>  V»y+±  =  1- 


Similarly,  if  =  - 


s}  +  j’ 


Sy+i  —  'S;+2.  =>  fpj+l  ~  1  =>  V’y+i.  —  1- 


It  also  follows,  using  (6.1),  that  the  inequalities  in  (ii)  and  (iii)  of  the  remarks  of  Section  V 
are  satisfied.  Also,  for  these  two  cases,  setting  Ty+  i.  =  1  —  0;  )_  ^  guarantees  that  Ty+  i.  =  0. 


Ty.fi  =  1  -  max(V»y,  0y+i). 
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The  flux  limiters  t/>  and  F  can  be  connected  to  those  used  in  upwind  schemes.  Sweby  15 
considers  an  upwind  Lax-Wendroff  scheme.  In  particular,  for  the  onc-dimcnsiona!  wave 


equation 


ut  +  aux  —  0, 


the  numerical  solution  is  obtained  with 


where 


Uy+1  —  Uj  —  —  A. 


v  —  a——,  0  <  <p(r)  <  2, 

Ax 


^(1  A 


<Pj(UJ  +  l  ~  Uj) 


and  the  backward  difference  operator  A_  is  defined  by 

A  __  uy  =  Uj  -  Uj  i . 

If  (6.2)  is  rewritten  as  a  central-difference  scheme,  then 

Uy  +  1  -  Uj  -  j(uj+ 1  -  Uj  i )  +  ^-[(1  -  <pj)  (uj+ j  -  u;)  -  (1  -  '-Pj-i)  [u}  -  u,  ,); 

V  2 

+  2  [^y(“/+i  -  Uj)  -  y?y-i  (uy  -  Uj_i))  - 

(< 

By  dropping  the  i/2  term  in  (6.3)  and  changing  to  a  semidiscrete  formulation,  we  get 
d  u  (  ^ 

+  ^  [(1  -  <Pj){uj+ 1  -  Uj)  -  (1  -  ¥>y_l)(uy  -  Uj-!]  . 


Pi  =  1  -  V'y+fc. 

then  the  second-difference  dissipation  term  has  the  same  form  as  the  one  presented  in 
Section  V. 

To  complete  the  connection  between  the  limiters  for  the  central-difference  and  upwind 
schemes,  we  compare  their  behavior.  For  the  central-difference  case, 

0<V'<1,  0  <  <p  <  1, 


whereas  for  the  upwind  case,  0  <  <p  <  2.  Furthermore,  for  the  central-difference  limiter, 

^(r)=V?(r) 
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means  that  it  does  not  matter  if  Au^i  >  Ait;_i  or  AuJ+i  <  Au;_i  (i.e.,  the  sign  of  a 
does  not  affect  the  scheme),  as  opposed  to  the  upwind  limiter  where 


In  Figure  2  the  Sweby  diagram  for  the  upwind  van  Leer  flux  limiter  and  the  central- 
difference  van  Leer  flux  limiter  is  shown.  For  r  <  1,  the  limiters  are  the  same.  For  r  >  1, 
the  upwind  version  continues  until  <p  =  2,  while  the  central-difference  version  returns  to 
zero.  In  Sweby’s  paper  the  flux  limiters  are  not  allowed  to  decrease  when  r  >  1.  However, 
there  is  no  difficulty  from  such  behavior.  In  both  cases  <p(r)  =  0  when  r  <  0  so  that  the 
limiter  is  turned  ofT  for  extrema.  We  note  that  other  limiters  besides  that  of  van  Leer  can 
be  used  to  get  switches  for  central-difference  schemes. 

For  systems  of  equations  we  use  a  scalar  limiter.  Using  the  matrix  form  of  the  dissi¬ 
pation  it  is  easy  to  implement  different  limiters  for  different  characteristic  variables. 

In  terms  of  the  pressure  the  switch  becomes 


/  1  “  r 

*“rr m 


Apj-k 

r  =  — - — ,  a  >  0 


Vb  = 


_  \P:+i  -Pj\~  (Pj  -  Pj-i)  sgn  ( P j 4  i  -  Pj) 


\Pj+ 1  ~Pj\  +  I  Pj  -  Pj- 1|  +  e 


a  >  0 


0j+i.  =  max(Vb,  ipj+i), 

and  e  =  0(Ax2)  to  prevent  a  zero  denominator  for  constant  pressure  regions.  Consider 
also  the  wave  speed  a  <  0.  Then, 


Vb  = 


P}  + 1  2py  +  Pj~  i 


|P;  +  i  Pj\  +  \Pj  Pj  i|  +  e  [  sgn(py_i  —  p;),  a<0 


sgn(p/+i  -  Pj),  a  >  0 


We  use  a  conservative  approach  and  take 


jPj  +  t  ~  2P;  +  Pj- 1! 

|P>+i  -  Pj  I  +  I  Pj  -  Pj  —  i  I  +  c ' 


Notice  that  this  switch  is  very  similar  to  (2.9)  for  the  original  dissipation  model  of  the 
Runge-Kutta  scheme.  There  is  only  a  minor  change  in  the  denominator.  However,  with 
this  change  and  the  factor  1/2  in  front  of  the  second-difference  dissipation  term,  the  scalar 


equation  becomes  first-order  upwind  near  shocks.  In  the  case  of  the  original  xpj  we  find 
that  xjjj  ~  .05  near  shock  waves  in  transonic  flows.  One  may  require  different  parameters 
for  the  Runge-Kutta  scheme  to  ensure  stability. 

We  now  no  longer  have  a  free  parameter  for  the  second-difference  dissipation.  We  also 
usually  use  F  =  1  —  20  so  that  the  fourth-difference  dissipation  is  cut  off  for  0  >  1/2.  The 
only  free  parameter  is  the  coefficient  of  the  fourth-difference  term. 

VII.  Numerical  Algorithm 

The  majority  of  the  numerical  results  presented  in  this  paper  were  obtained  with  a 
Navier-Stokes  code  developed  by  the  authors,  which  is  based  on  the  explicit  multistage 
time-stepping  schemes  of  [3]  and  [16].  This  class  of  schemes  is  currently  in  widespread 
use  for  solving  the  Euler  equations.  In  references  [17-20],  these  schemes  were  extended 
to  allow  solution  of  the  compressible  Navier-Stokes  equations.  Significant  improvements 
in  numerical  efficiency  were  introduced  in  [9],  [10],  and  [21  ] .  In  the  code  of  Swanson  and 
Turkel  a  cell-centered,  finite-volume  method  is  employed  to  obtain  centered  type  difference 
approximations  for  the  flow  equations.  Such  a  method  provides  flexibility  in  treating 
arbitrary  geometries  and  different  grid  topologies,  since  no  special  treatment  is  required  in 
the  vicinity  of  singular  points  or  lines.  The  scheme  is  second-order  accurate  in  space  for 
sufficiently  smooth  meshes  (see  [21-22]  for  definition  of  sufficiently  smooth).  A  modified 
five  stage  Runge-Kutta  scheme  is  used  for  the  time  integration  to  obtain  a  steady-state 
solution.  There  is  a  weighted  evaluation  of  the  artificial  dissipation  terms  on  the  first, 
third,  and  fifth  stages,  which  provides  a  large  parabolic  stability  limit.  The  physical 
viscous  terms  are  computed  on  the  first  stage  only  and  frozen  for  the  remaining  ones;  this 
does  not  compromise  the  stability  characteristics  of  the  scheme.  The  spatial  and  temporal 
differencing  are  decoupled.  Thus,  the  numerical  algorithm  is  independent  of  time  step 
and  amenable  to  steady-state  convergence  acceleration  techniques.  These  methods  include 
local  time  stepping  (a  preconditioning  for  the  system  of  difference  equations),  variable 
coefficient  implicit  residual  smoothing,  and  multigrid.  Implicit  residual  smoothing  is  just 
a  mathematical  step,  applied  at  each  stage  of  the  explicit  scheme,  to  extend  the  local 
stability  range.  The  multigrid  method  involves  cycling  through  a  sequence  of  successively 
coarser  grids  and  relying  upon  effective  high  frequency  damping  for  rapid  removal  of  the 
errors  in  the  fine  grid  solution. 


17 


VIII.  Applications 


Two  airfoil  flow  problems  are  considered  here  to  demonstrate  the  benefits  of  using  a 
central-difference  scheme  with  a  matrix  numerical  dissipation  model.  The  first  problem 
concerns  inviscid  flow  over  an  NACA  0012  airfoil.  At  a  free-stream  Mach  number  (Moo)  of 
0.8  and  angle  of  attack  (a)  of  1.25  degrees,  a  fairly  strong  transonic  shock  wave  occurs  on 
the  upper  surface,  while  a  much  weaker  one  appears  on  the  lower  surface.  On  representative 
inviscid  meshes,  schemes  based  on  central  differencing  capture  the  upper  surface  shock 
reasonably  well,  but  they  smear  significantly  the  lower  surface  shock.  The  second  problem 
involves  transonic  turbulent  flow  over  a  RAE  2822  airfoil.  In  this  case  the  free-stream 
Mach  number  is  0.73,  the  angle  of  attack  is  2.79  degrees,  and  the  Reynolds  number  based 
on  chord  (Reoo)  is  6.5  x  106.  For  such  transonic  viscous  flows  small  differences  in  the  shock 
strength  can  result  in  noticeable  changes  in  the  lift  and  drag  of  the  airfoil.  Thus,  both 
of  these  problems  can  provide  a  reasonable  measure  of  the  performance  of  the  artificial 
dissipation  model,  especially  near  shock  waves. 

A  C-type  mesh  consisting  of  224  cells  around  the  airfoil  (160  cells  on  the  airfoil) 
and  32  cells  normal  to  the  airfoil  was  used  for  the  first  problem.  The  outer  boundary 
of  the  finite  domain  was  placed  20  chords  away  from  the  airfoil,  so  as  to  not  produce 
significant  effects  on  the  solutions  (see  [23]).  The  normal  mesh  spacing  at  the  surface  was 
approximately  0.01  chords.  The  mesh  was  clustered  at  the  leading  and  trailing  edges  of  the 
airfoil  in  order  to  improve  resolution  of  the  steep  gradients  occurring  in  these  regions.  In 
Figure  3  the  surface  pressure  distributions  computed  with  the  scalar  and  matrix  dissipation 
models  are  shown.  There  is  a  discernible  improvement  in  the  sharpness  of  the  shock 
waves  using  the  matrix  model,  especially  for  the  one  on  the  lower  airfoil  surface.  The 
lift  and  drag  coefficients  obtained  using  the  two  models,  along  with  the  values  from  a 
high  density  mesh  calculation,  are  compared  in  Table  I.  The  values  associated  with  the 
matrix  model  nearly  match  those  corresponding  to  the  fine  mesh  solution.  The  convergence 
histories  in  terms  of  number  of  multigrid  cycles  are  also  displayed  in  Figure  3.  Convergence 
is  measured  by  the  logarithm  of  the  root  mean  square  of  the  residual  of  the  continuity 
equation.  Th*>  convprgence  rates  obtained  with  the  scalar  ^nd  matrix  models  are  0.819 
and  0.888,  respectively.  The  matrix  model  result  required  about  55  epu  sec  on  a  Cray  II 
computer  for  125  multigrid  cycles  with  the  solution  grid  This  processing  time  represents 
a  15  percent  increase  in  the  time  needed  compared  with  the  scalar  dissipation  model.  It 
may  be  possible  to  improve  the  convergence  rate  with  the  matrix  model  by  altering  the 
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coefficients  of  the  time-stepping  scheme  and/or  the  type  of  implicit  residual  smoothing 
operator  employed. 

The  following  set  of  C-type  meshes  were  used  in  solving  the  second  problem  with 
the  two  dissipation  models:  l)  160  x  32  with  128  cells  on  the  airfoil,  2)  320  x  64  with 
256  cells  on  the  airfoil,  3)  640  x  128  with  512  cells  on  the  airfoil.  Each  successively  coarser 
mesh  was  generated  by  eliminating  every  other  mesh  line  in  both  coordinate  directions  of 
the  finer  mesh.  Again  the  outer  boundary  of  the  domain  was  located  20  chords  from  the 
airfoil.  The  normal  spacing  at  the  surface  for  the  finest  mesh  was  7.5  x  10~6  chords.  In 
Figure  4  a  comparison  is  made  between  the  experimental  data  of  [24]  and  the  predicted 
distributions  of  pressure,  upper  surface  skin  friction,  and  upper  surface  boundary-layer 
displacement  thickness  obtained  on  the  160  x  32  mesh.  The  skin-friction  coefficient  is  the 
wall  shear  stress  nondimensionalized  by  the  dynamic  pressure  at  the  edge  of  the  viscous 
layer.  Velocity  profiles  (scaled  by  the  boundary-layer  edge  velocity  ue)  at  two  stations  on 
the  upper  airfoil  surface  are  also  presented  in  Figure  4.  There  is  marked  improvement  in 
the  results  obtained  with  the  matrix  model,  as  verified  by  their  reasonable  agreement  with 
the  data.  The  convergence  histories  given  in  Figure  4  indicate  that  the  final  residual  level 
is  somewhat  higher  with  the  matrix  model.  However,  the  terminal  rate  of  convergence  is 
nearly  the  same  with  both  dissipation  models.  Finally,  Figures  5  and  6  show  the  pressure 
and  skin-friction  distributions  computed  on  the  coarse  mesh  with  second-order  upwind  and 
“third-order”  upwind  biased  forms  of  Roe’s  scheme  [2] ,  which  has  been  shown  to  have  low 
levels  of  dissipation.  The  pressures  obtained  with  the  “third-order”  form  are  very  close  to 
those  calculated  with  the  matrix  model.  The  upwind  skin-friction  solutions  exhibit  slightly 
better  agreement  with  the  experimental  data  upstream  of  the  shock. 

Numerical  results  for  the  two  finer  meshes  are  presented  in  Figures  7  and  8.  For  the 
320  x  64  mesh,  a  slightly  stronger  shock  is  predicted  using  the  matrix  dissipation  model. 
Otherwise,  the  solutions  determined  with  the  different  dissipation  models  are  nearly  the 
same.  The  character  of  the  convergence  behavior  with  the  two  models  is  about  the  same 
as  for  the  coarse  mesh.  On  the  finest  mesh,  the  application  of  the  models  gives  essentially 
the  same  results.  In  Figure  9  the  variation  of  the  computed  lift  and  total  drag  coefficients 
with  the  reciprocal  of  the  number  of  mesh  cells  is  plotted  for  each  dissipation  model. 
The  curves  corresponding  to  the  matrix  model  are  nearly  linear.  A  linear  curve  indicates 
second-order  accuracy  for  the  full  range  of  meshes  being  considered.  The  numerical  values 
for  the  components  of  the  drag  coefficient  (form  and  friction  contributions)  along  with  the 
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lift  coefficient  are  given  for  all  the  viscous  computations  in  Table  II.  Notice  that  the  values 
obtained  with  the  matrix  model  on  the  320  x  64  grid  are  very  close  to  those  calculated 
with  the  scalar  model  on  the  640  x  128  grid.  Such  a  reduction  in  the  mesh  required  for 
acceptable  accuracy  is  highly  desirable,  especially  for  three-dimensional  flow  problems. 
Some  applications  to  three-dimensional  viscous  flows  are  presented  in  Turkel  and  Vatsa 
[25],  It  should  be  emphasized  that  accurate  prediction  of  lift  and  drag  is  very  important 
in  the  design  of  aircraft,  and  thus,  a  good  estimate  of  these  quantities  for  an  infinitely  fine 
mesh  is  needed. 

In  [26]  it  is  shown  that  the  TVD  switch  (6.4)  allows  converged  solutions  for  hypersonic 
flow  where  the  standard  switch  (2.9)  does  not. 

IX.  Concluding  Remarks 

We  have  thus  shown  that  the  second-difference  artificial  dissipation  is  equivalent  to 
using  a  flux  limiter,  and  hence,  central-difference  schemes  are  not  any  more  “artificial” 
than  upwind  schemes.  The  central-difference  scheme  is  slightly  more  dissipative  for  two 
reasons.  First,  =  max(V>y ,  V’y+i)  while  for  an  upwind  scheme  V'y+i  is  equal  to  either 

x/jj  or  0y+i,  depending  on  the  direction  of  the  wind.  Second,  we  insist  that  xpj  be  positive 
(i.e.,  <py  <  1)  while  upwind  limiters  allow  <pj  >  1  (i.e.,  in  some  cases  we  can  have  negative 
viscosity  but  still  be  TVD).  However,  to  compensate  for  this  slight  increase  in  dissipation 
central-difference  schemes  are  simpler  to  program  and  require  less  computer  time  per  time 
step,  and  work  well  with  multigrid  acceleration  techniques. 

In  addition  the  central-difference  schemes  have  a  free  parameter  in  conjunction  with 
the  fourth-difference  dissipation.  This  dissipation  is  needed  to  approach  a  steady  state 
and  has  nothing  to  do  with  TVD  properties.  In  fact  the  fourth-difference  contribution  is 
set  equal  to  zero  near  local  extrema.  For  time-dependent  flows  one  can  set  this  dissipation 
identically  to  zero.  On  the  other  hand  TVD  properties  do  not  necessarily  imply  a  rapid 
convergence  to  the  steady  state. 

In  summary,  a  formulation  for  a  numerical  dissipation  model  that  makes  a  central- 
difference  scheme  closely  resemble  an  upwind  scheme  near  flow  discontinuities  has  been 
described.  A  theorem  has  been  proven  that  gives  sufficient  requirements  for  this  type  of 
dissipation  model  to  satisfy  the  TVD  property  for  a  scalar  equation.  Flux  limiter  functions 
have  been  presented  for  this  form  of  dissipation  model.  For  a  system  of  equations  a  matrix¬ 
valued  dissipation  is  introduced.  Solutions  of  the  Euler  and  Navier-Stokes  equations  for 
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airfoil  flows  have  been  obtained  using  the  matrix  dissipation  model.  These  results  have 
demonstrated  noticeable  improvements  in  accuracy  in  smooth  regions  of  the  flow  field 
as  well  as  near  shock  waves.  There  is  a  15  percent  increase  in  computational  time  for 
explicit  multistage  schemes  when  this  model  is  used.  However,  the  improved  accuracy  has 
permitted  a  significant  reduction  in  the  number  of  mesh  points  required.  Such  behavior  of 
a  scheme  can  have  a  dramatic  effect  on  the  necessary  mesh  size  for  three-dimensional  flow 
calculations. 

Finally,  it  is  important  to  emphasize  the  different  principal  objectives  associated  with 
matrix-valued  dissipation  and  the  TVD  switch.  The  purpose  of  the  matrix  form  of  the 
numerical  dissipation  model  is  to  apply  the  appropriate  scaling  of  the  dissipation  in  each 
flow  equation,  yielding  a  reduction  in  the  amount  of  dissipation  being  introduced  and 
improved  accuracy.  The  TVD  switch  plays  a  somewhat  opposite  role  in  that  it  causes 
more  dissipation  to  be  added  in  order  to  prevent  overshoots,  and  thus,  allows  convergence 
and  provides  robustness  in  solving  high  speed  flow  problems.  The  combination  of  the  two 
should  give  more  dissipation  near  shocks  and  less  dissipation  in  smooth  regions,  hence 
giving  better  accuracy  in  all  regions. 
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Table  I 

Lift  and  drag  coefficients  for  NACA  0012  airfoil, 
Moo  =  0.8,  a  =  1.25° 


Case 

ci 

Cd 

Scalar  dissipation  model 
(224  x  32  mesh) 

0.3628 

0.0231 

Matrix  dissipation  model 
(224  x  32  mesh) 

0.3591 

0.0227 

Scalar  dissipation  model 
(640  x  64  mesh) 

0.3577 

0.0228 

C|  -  lift  coefficient 
Cd  -  drag  coefficient 


i 

« 


i 
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Table  II 

Lift  and  drag  coefficients  for  RAE  2822  airfoil, 
Moo  =  0.73,  a  =  2.79°,  Re^  =  6.5  x  106 


Scalar  dissipation  model 


Mesh 

Cl 

c<tp 

cdf 

Cdtot 

160  x  32 

0.8081 

0.0128 

0.0045 

0.0173 

320  x  64 

0.8395 

0.0120 

0.0054 

0.0174 

640  x  128 

0.8544 

0.0123 

0.0055 

0.0178 

Matrix  dissipation  model 


160  x  32  0.8296  0.0124  0.0050  0.0174 

320  x  64  0.8514  0.0123  0.0055  0.0178 

640  x  128  0.8588  0.0124  0.0055  0.0179 

Second-order  upwind 


160  x  32  0.8220  0.0124  0.0054  0.0178 


ci  -  lift  coefficient 
Cdp  -  pressure  drag  coefficient 
Cdf  -  friction  drag  coefficient 
cd<ot  “  total  drag  coefficient 


Figure  2  Sweby  diagram 


Cycles 


(  a  ) 


(  b  ) 


Figure  3  Central-difference  scheme  results  for  inviscid  flow  over  NACA  0012  airfoil 
(Moo  =  0.80,  a  =  1.25  degrees): 

(a)  Surface  pressure  coefficient,  (b)  Convergence  history 
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(  a  )  (  b  ) 


Figure  4  Central-difference  scheme  results  for  turbulent  flow  over  RAE  2822  airfoil 
(160  x  32  mesh,  =  0.73,  a  =  2.79  degrees,  Re^  =  6.5  x  106): 

(a)  Surface  pressure  coefficient,  (b)  Convergence  history 


(  C  )  (  d  ) 


Figure  4c  -  4d  Upper  surface  skin-friction  coefficient  (cf)  and  boundary-layer  displacement 
thickness  (160  x  32  mesh) 


30 


50  x  10'3 


(  e  ) 


(  f  ) 


Figure  4e  -  4f  Velocity  profiles  (160  x  32  mesh) 


7  x  10-3 


(  a  ) 


(  b  ) 


Figure  5  Second-order  Roe  scheme  results  for  turbulent  flow  over  RAE  2822  airfoil 
(160  x  32  mesh,  =  0.73,  a  =  2.79  degrees,  Re^  =  6.5  x  106): 

(a)  Surface  pressure  coefficient,  (b)  Upper  surface  skin-friction  coefficient 
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Figure  6  “Third-order”  Roe  scheme  results  for  turbulent  flow  over  RAE  2822  airfoil 
(160  x  32  mesh,  =  0.73,  a  =  2.79  degrees,  Re^  =  6.5  x  106): 

(a)  Surface  pressure  coefficient,  (b)  Upper  surface  skin-friction  coefficient 


Cycles 


(  a  ) 


(  b  ) 


Figure  7  Central-difference  scheme  results  for  turbulent  flow  over  RAE  2822  airfoil 
(320  x  64  mesh,  =  0.73,  a  -  2.79  degrees,  Re^  =  6.5  x  106): 

(a)  Surface  pressure  coefficient,  (b)  Convergence  history 
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Figure  7c  -  7d  Upper  surface  skin-friction  coefficient  (cf)  and  boundary-layer  displacement 
thickness  (320  x  64  mesh) 
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(f  ) 

Figure  7e  -  7f  Velocity  profiles  (320  x  64  mesh) 

Figure  8  Central-difference  scheme  results  for  turbulent  flow  over  RAE  2822  airfoil 
(640  x  128  mesh,  =  0.73,  a  =  2.79  degrees,  Re^  =  6.5  x  106): 

(a)  Surface  pressure  coefficient,  (b)  Convergence  history 
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Figure  8c 


7  x  10'3 


Upper  surface  skin-friction  coefficient  (6-10  x  128  mesh) 
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1/N  x  10~4 


Figure  9a  Variation  of  lift  coefficient  with  reciprocal  of  number  of  mesh  points 
(RAE  2822  airfoil,  =  0.73,  a  -  2.79  degrees,  Re^  6.5  x  10r,l 
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